#!/usr/bin/python
# coding: utf-8

import sys
import logging
import subprocess
from math import sqrt, isnan
from itertools import izip, groupby, chain
from itertools import permutations, product
from math import log, pi, sqrt, copysign
from random import shuffle

try:
    from scipy.integrate import quad
    from scipy.special import k0
    from scipy.stats import norm
    
    NORM_DIST = norm(0.0, 1.0)
except ImportError:
    pass

from optparse import Option, OptionParser

from bessel import BesselNormValues

ALPHABET = 'actg'
PERIOD   = 3

def inv_xrange(n):
    for i in xrange(n):
        yield n - i - 1

MAX_BESSEL_VALUE = 500.0

def AbsBessel2Norm(abs_x):
    if abs_x <= 4:
        integr_value, err = quad(k0, 0, abs_x)
        prob = 0.5 - integr_value / pi
    else:
        if abs_x > MAX_BESSEL_VALUE:
            logging.warn('Bessel2Norm: abs(x) = %s', abs_x)
            abs_x = MAX_BESSEL_VALUE
        integr_value, err = quad(k0, abs_x, MAX_BESSEL_VALUE*2)
        prob = 2 / pi * integr_value
    assert err < 0.001
    return NORM_DIST.isf(prob)

def SignedSafeCall(func, x):
    if x != 0:
        res = func(abs(x))
        if x < 0 and res > 0:
            return -res
        else:
            return res
    else:
        return 0.0

def generate_Bessel2Norm():
    d = {}
    for n in xrange(1, int(MAX_BESSEL_VALUE)):
        precision = 10**(2 - int(log(n)/log(10)))
        assert precision > 0
        for p in xrange(0, precision):
            x = n - p / float(precision)
            print x
            d[x] = AbsBessel2Norm(x)
    return d

Bessel2Norm = lambda x: SignedSafeCall(AbsBessel2Norm, x)

def Bessel2NormFix(x):
    ai = min(1000, abs(int(x*10)));
    if (x >= 0):
        return BesselNormValues[ai]
    else:
        return -BesselNormValues[ai]

def clog(f):
    if f > 0:
        return f * log(f)
    else:
        return 0.0

def mutual_inf_mx(m1, m2):
    I = 0
    for j in xrange(PERIOD):
        y1, y2 = 0, 0
        for c in ALPHABET:
            w1 = m1[j][c]
            w2 = m2[j][c]
            y1 += w1
            y2 += w2
            I += clog(w1) + clog(w2) - clog(w1 + w2)
        I += clog(y1 + y2) - clog(y1) - clog(y2)
    return I

def inf_mx(mx):
    I = 0
    y = dict((j,sum(mx[j][c] for c in ALPHABET)) for j in xrange(PERIOD))
    x = dict((c,sum(mx[j][c] for j in xrange(PERIOD))) for c in ALPHABET)
    L = float(sum(x.itervalues()))
    assert L == sum(y.itervalues())    
    for j in xrange(PERIOD):
        for c in ALPHABET:
            I += clog(mx[j][c])
    for j in xrange(PERIOD):
        I -= clog(y[j])
    for c in ALPHABET:
        I -= clog(x[c])
    I += clog(L)
    return I
    
def normalize_mx(mx):
    y = dict((j,sum(mx[j][c] for c in ALPHABET)) for j in xrange(PERIOD))
    x = dict((c,sum(mx[j][c] for j in xrange(PERIOD))) for c in ALPHABET)
    L = float(sum(x.itervalues()))
    assert L == sum(y.itervalues())    
    nmx = {}
    for j in xrange(PERIOD):
        nmx[j] = {}
        for c in ALPHABET:
            p = y[j] * x[c] / L / L
            nmx[j][c] = p and (mx[j][c] - p * L) / sqrt(L * p * (1 - p)) or 0.0
    return nmx

def norm_dist_mx(nmx1, nmx2):
    D = 0.0
    for j in xrange(PERIOD):
        for c in ALPHABET:
            D += Bessel2NormFix(nmx1[j][c] * nmx2[j][c])
    return D

def create_mx(seq):
    d = dict((j, dict((c,0) for c in ALPHABET)) for j in xrange(PERIOD))
    for i, c in enumerate(seq):
        if c in ALPHABET:
            d[i%PERIOD][c] += 1
    #print seq
    return d

def create_norm_mx(seq):
    return normalize_mx(create_mx(seq))

def print_mx(mx):
    for c in ALPHABET:
        print "\t".join("%.3f" % mx[j][c] for j in xrange(PERIOD))


def print_dp_mx(dp_mx):
    TERM_WIDTH = 180
    COL_WIDTH  = 12
    COL_PER_LINE = TERM_WIDTH / COL_WIDTH
    N = len(dp_mx[0])
    L = len(dp_mx)
    for b in xrange(0, L, COL_PER_LINE):
        for i in xrange(N):
            print " ".join("%7.2f %i" % dp_mx[j][i] for j in xrange(b, min(L,b+COL_PER_LINE)))
        print


def create_dp_mx(options, base_mx, seq):
    reward     = options.reward
    penalty    = options.penalty
    win_length = options.win_length
    length     = len(seq)
    seq        = seq.lower()
    
    bmxs = [dict(((i-j)%PERIOD, n) for i,n in base_mx.iteritems()) for j in xrange(PERIOD)]
    #print bmxs
    #for j in xrange(PERIOD):
    #    print '---------'
    #    print_mx(bmxs[j])
    
    dp_mx = [[(0,i) for i in xrange(PERIOD)]]
    #print_dp_mx(dp_mx)
    
    for x in xrange(0, length-win_length-1):
        #print x, "################################"
        mxs = [create_norm_mx(seq[x+j:x+j+win_length]) for j in xrange(PERIOD)]
        #print_mx(bmxs[x % PERIOD])
        #print '***'
        #print_mx(mxs[0])
        #print '***'
        deltas = [norm_dist_mx(bmxs[x%PERIOD], mxs[j]) for j in xrange(PERIOD)]
        #print 'deltas: ', ' | '.join(map(str, deltas))
        new_col = [
            max(
                (
                    (
                        max(dp_mx[x][j][0] + (i==j and deltas[j]-reward or -penalty), 0),
                        j
                    )
                    for j in (k%PERIOD for k in xrange(i, i+PERIOD))
                ),
                key = lambda x:x[0],
            )
            for i in xrange(PERIOD)
        ]
        dp_mx.append(new_col)
        #print "%5i: %s" % (x, ''.join("(%7.2f, %i), " % _i for _i in new_col))
        #print rewards
        #print_dp_mx(dp_mx)
    
    return dp_mx

def create_backtrace(options, dp_mx):
    max_reward = float('-inf')
    for i in xrange(len(dp_mx)):
        for j in xrange(PERIOD):
            reward = dp_mx[i][j][0]
            if reward > max_reward:
                max_reward = reward
                max_i, max_j = i, j
    
    #print 'RES:', max_reward, max_i, max_j

    backtrace = []
    start_pos = 0
    if max_reward > 0:
        shift = max_j
        for i in inv_xrange(max_i+1):
            reward, shift = dp_mx[i][shift]
            #print i, reward, shift
            backtrace.append(shift)
            if reward <= 0:
                break
        backtrace.reverse()
        start_pos = max_i - len(backtrace) + 1
        last_shift = backtrace[-1]
        backtrace.extend([last_shift] * ( ((len(backtrace) + options.win_length)/PERIOD)*PERIOD - len(backtrace) ) )
    return max_reward, start_pos, backtrace


def print_res(subseq, backtrace):
    if backtrace:
        last_shift = backtrace[0]
        i = 0
        for j, shift in enumerate(backtrace):
            diff = (PERIOD - abs(shift - last_shift)) % PERIOD
            for k in xrange(diff):
                i += 1
                sys.stdout.write('-')
                if i % PERIOD == 0:
                    sys.stdout.write('|')
            sys.stdout.write(subseq[j])
            i += 1
            if i % PERIOD == 0:
                sys.stdout.write('|')
            last_shift = shift
    sys.stdout.write('\n')


def sort_cycle(l):
    minindex = l.index(min(l))
    return l[minindex:]+l[:minindex]

def uniq(iterable, key=None):
    for key, grp in groupby(iterable, key=key):
        yield grp.next()

def iter_base_period_seq():
    perms = ((perm, sort_cycle(perm)) for perm in permutations(ALPHABET, PERIOD))
    for p, sp in uniq(sorted(perms, key=lambda i:i[1]), key=lambda i:i[1]):
        yield ''.join(p)
    
def iter_base_mx(win_length, parts_num=None):
    """
    Генерирует всевозможные варианты последовательностей длины win_length
    для создания базовых матриц с заданным уровнем точности:
      * iter_base_period_seq выдаёт все варианты последовательностей
        длины PERIOD из символов ALPHABET с учётом циклических перестановок
      * при parts_num = 1  - минимальная точность,
        при parts_num = len(ALPHABET) - средняя,
        при parts_num = win_length / PERIOD - максимальная (все варианты)
    """
    parts_num = parts_num or len(ALPHABET)
    assert win_length % PERIOD == 0
    if parts_num > win_length / PERIOD: parts_num = win_length / PERIOD
    base_seq_parts = set((tuple(sorted(p)) for p in product(iter_base_period_seq(), repeat=parts_num)))

    min_win_length_part = win_length / (PERIOD * parts_num) * PERIOD
    max_win_length_part = win_length - min_win_length_part * (parts_num - 1)
    win_length_parts = [min_win_length_part] * (parts_num - 1) + [max_win_length_part]

    for base_parts in base_seq_parts:
        seq = ""
        for base_seq, win_length_part in izip(base_parts, win_length_parts):
            seq += base_seq * (win_length_part/PERIOD)
        #print seq
        yield seq


def find_random(options, (seq_id, seq), parts_num=4):
    max_reward = 0
    best_cnt_shifts = 0
    best_mx = None
    for base_seq in iter_base_mx(options.win_length, parts_num):
        #print base_seq
        #print seq
        base_mx = create_norm_mx(base_seq)
        #print_mx(base_mx)
        #print 
        dp_mx = create_dp_mx(options, base_mx, seq)
        #print_dp_mx(dp_mx)
        reward, start_pos, backtrace = create_backtrace(options, dp_mx)
        cnt_shifts = max(0, len(list(uniq(backtrace))) - 1)
        print "REWARD =", cnt_shifts, "[%i:%i]\t"%(start_pos, start_pos+len(backtrace)), reward
        #subseq = seq[start_pos:start_pos+len(backtrace)]
        #print len(backtrace), backtrace
        #print subseq
        #print_res(subseq, backtrace)
        if reward > max_reward:
            max_reward = reward
            best_mx = base_mx
            best_cnt_shifts = cnt_shifts
    print "%s\t%s\t%s\t%s" % (seq_id, len(seq), best_cnt_shifts, max_reward)
    #print_mx(best_mx)

def run_cpp_finder(opts, base_seq, seq):
    p = subprocess.Popen(
        args = ["../cpp/shift_finder_dp -l%s -r%s -p%s -m%s -s%s" % (
                opts.win_length, opts.reward, opts.penalty, base_seq, seq,
            )
        ],
        stdout = subprocess.PIPE,
        shell = True,
    )
    p.wait()
    for line in p.stdout:
        pass
    return float(line.strip())
    

def find_random_ext(options, seq, parts_num=4):
    max_reward = None
    best_mx = None
    for base_seq in iter_base_mx(options.win_length, parts_num):
        #print base_seq
        #print seq
        reward = run_cpp_finder(options, base_seq, seq)
        #print "REWARD =", reward
        if reward > max_reward:
            max_reward = reward
            best_mx = create_norm_mx(base_seq)
    print "MAX_REWARD =", max_reward
    #print_mx(best_mx)
    return max_reward

def pshuffle(seq):
    assert len(seq) % PERIOD == 0
    indices = [[seq[i] for i in xrange(j % PERIOD, len(seq), PERIOD)] for j in xrange(PERIOD)]
    for l in indices:
        shuffle(l)
    return ''.join(chain(*izip(*indices)))

def find_mk(options, seq, iter_num=100):
    base_mx    = create_norm_mx('act'*(options.win_length/PERIOD))
    print_mx(base_mx)
    print norm_dist_mx(base_mx, base_mx)
    dp_mx = create_dp_mx(options, base_mx, seq)
    max_reward, start_pos, backtrace = create_backtrace(options, dp_mx)
    subseq = seq[start_pos:start_pos+len(backtrace)]
    print subseq

    dp_mx = create_dp_mx(options, base_mx, subseq)
    max_reward, start_pos, backtrace = create_backtrace(options, dp_mx)

    print list(backtrace)
    print list(uniq(backtrace))
    print_dp_mx(dp_mx)
    print "REWARD =", max_reward
    #print
    rewards = []
    for i in xrange(iter_num):
        rseq = pshuffle(subseq)
        rdp_mx = create_dp_mx(options, base_mx, rseq)
        reward, start_pos, backtrace = create_backtrace(options, rdp_mx)
        print list(uniq(backtrace))
        #print_dp_mx(rdp_mx)
        print "REWARD =", reward
        rewards.append( reward )
    
    M = sum(rewards) / iter_num
    D = sum((reward - M)**2 for reward in rewards ) / iter_num 
    Z = (max_reward - M) / sqrt(D)
    print Z, M, D

def find_mk_ext(options, seq, iter_num=100):
    raise Exception('Not impemented!')


def find_strict(options, seq):
    base_mx    = create_norm_mx('act'*(options.win_length/PERIOD))
    #print_mx(base_mx)
    dp_mx = create_dp_mx(options, base_mx, seq)
    #print_dp_mx(dp_mx)
    max_reward, start_pos, backtrace = create_backtrace(options, dp_mx)
    print max_reward
    print start_pos, len(backtrace), len(seq)
    subseq = seq[start_pos:start_pos+len(backtrace)]
    #print len(backtrace), backtrace
    #print subseq
    print_res(subseq, backtrace)


if __name__ == '__main__':
    parser = OptionParser(usage='%prog [-r REWARD -p PENALTY -l WIN_LENGTH -s SEQUENCE]')
    parser.add_option(Option(
        '-r', '--reward',
        dest = 'reward',
        type = int,
        default = 1,
    ))
    parser.add_option(Option(
        '-p', '--penalty',
        dest = 'penalty',
        type = int,
        default = 3,
    ))
    parser.add_option(Option(
        '-l', '--window-length',
        dest = 'win_length',
        type = int,
        default = 36,
    ))
    parser.add_option(Option(
        '-s', '--sequence',
        dest = 'sequence',
    ))
    opts, args = parser.parse_args()
    if opts.win_length % PERIOD != 0:
        print >> sys.stderr, "WIN_LENGTH is not divided by %s" % PERIOD
        sys.exit(1)
        
    if opts.sequence:
        #find_strict(opts, opts.sequence)
        #find_mk(opts, opts.sequence, 10)
        find_random_ext(opts, opts.sequence)
    else:
        for line in sys.stdin:
            id, line = line.strip().split()
            #find_strict(opts, line)
            #find_mk(opts, line, 100)
            find_random(opts, (id, line), parts_num=2)
            #find_random_ext(opts, line, parts_num=4)
            #find_mk_ext(opts, line, 10)
